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Properties of a sequence of positive operators defined by the Widder 
Laplace inversion formula are studied in order to obtain practical 
methods for the inversion of the Laplace transform, practical error 
formulae, and useful approximations to given functions. The approx- 
imation procedure retains essential structural characteristics of the 
original function, e.g., nonnegativity , monotonicity, and convexity. 
Thus a distribution function is approximated by distribution functions. 
Enhancement techniques are provided for the improvement of accuracy 
for a given order of approximation. The methods are illustrated by 
applications to renewal theory and to the covariance and recovery 
functions of telephone traffic theory. 

I. INTRODUCTION 

The Laplace transform occurs frequently in investigations of queueing 
theory and telephone traffic models in which it usually represents a 
probability distribution function. Although the mean and variance of 
the distribution can be readily obtained from the transform, there are 
many investigations in which the distribution itself is needed; in par- 
ticular, good analytic and numerical approximations for the comple- 
mentary distribution when the argument is large. This is the case, for 
example, when studying waiting times of queues, time delays of work 
through a computer system, and delays of message progress through data 
networks. 

Numerical methods which have been made thus far 16 - 18 concentrate 
on accurate numerical approximation on some interval [0, T], the diffi- 
culty of accurate inversion increasing with increasing T. Methods de- 
pending on Gauss-Legendre quadrature applied to the defining Laplace 
integral with subsequent interpolation are discussed in Ref. 19. These 
methods require the solution of large order linear systems whose matrices 
are severely ill-conditioned; thus they can bog down in meaningless 
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calculations. Much ingenuity has been used in specific cases to circum- 
vent this problem. Asymptotic formulae may sometimes be used to ap- 
proximate the complementary distribution for large argument; however, 
in many practical cases, good accuracy was obtained only when the 
argument was so large that the corresponding probabilities were too small 
to be of practical significance. One of the methods of this paper, namely, 
the a enhancement procedure, specifically attacks this problem by im- 
itating the exponential decay of the original in [T, <»] while simulta- 
neously providing accurate approximation in [0, T\. The transition region 
is sufficiently well approximated for most practical uses. 

The well-known Laplace inversion formula of Widder 1 - 3 has not been 
actively used in practical work. It has been the experience of the author 
that an investigation of the Widder formula qua functional transfor- 
mation can provide useful practical techniques for inversion and also 
inequalities and limit relations between the approximations and the 
original function. Accordingly, it is the object of this paper to study the 
properties of a sequence of positive operators defined by the Widder 
formula in order to obtain practical methods for Laplace inversion, 
practical error formulae, and useful approximations to given functions. 

In II the Widder inversion formula is obtained and a sequence of 
positive operators, L n , which form the subject of the paper, are intro- 
duced. The L n map a function f(t){t > 0) to a sequence of functions /„ (t) 
= L n f which converge uniformly on [0, «] to /(£). This viewpoint enables 
one to study the approximation characteristics of the sequence f n (t), thus 
providing a means of approximating a given f(t) besides effecting the 
inversion of its Laplace transform, /(s). Several representations are given 
for f n (t) in terms of f(t). 

In III properties of the sequence \f n (t)\ are developed which show that 
it possesses many desirable characteristics. In many applications it is 
preferable that the approximating functions globally imitate the original 
function in qualitative structural features rather than to the attainment 
of very high numerical accuracy. Thus if the original function lies be- 
tween zero and one, is monotone decreasing, and is convex, then these 
same properties would be desired in the approximation. It is shown that 
the approximating sequence, \f n (t)\, does retain those properties. A re- 
cursion relation for f n +i(t) in terms of f n (t), and a generating function 
for the sequence are also given, thus making the computation of higher 
approximations possibly more convenient than the direct application 
of the representation formulae themselves. A useful feature of the f n (t) 
is that, when f(t) is convex, they satisfy /„(£) > f(t). 

Part IV develops error bounds and pointwise error estimates. The 
results in terms of f(t) reflect the use of the technique for approximation; 
on the other hand, the pointwise estimate of error in terms of f n (t) is 
especially useful for the inversion problem since then f(t) is not available. 

670 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1978 



It is also shown that the successive approximations f (t), f x {t), f 2 (t), . . . 
are uniformly better for each t if f(t) > 0. 

In practical use the initial member of the sequence, fo(t), is not an 
accurate approximation to f(t) for t not in the neighborhoods of zero and 
infinity. Additionally the sequence \f n (t)\ does not converge rapidly in 
n. Consequently one must go far out in the sequence to obtain adequate 
accuracy. Part V treats this problem. A modification, /„,„(£), of f n (t) is 
introduced depending on a parameter a for which, by appropriate choice 
of a, fo,„(t) is a much improved approximation to f(t) than fo(t); the 
rapidity of convergence of the sequence \f n , a (t)\ is not improved over that 
of the unmodified sequence. However, it has been found that good ac- 
curacy is obtained by use of fo, a {t) or fi t „(t) as is demonstrated in the 
examples on covariance and recovery functions given in this paper. 

In many applications, especially to complementary distribution 
functions, the behavior of f(t) for large t must be accurately reproduced. 
The approximations /„,„(£) accomplish this especially when a is related 
to the decrement of an exponential majorant. For functions which are 
exponentially small at infinity, the f n (t) do not adequately reproduce 
the decay of f(t). 

Many of the desirable features of the original method are still retained 
by this modification. The concept of convexity with exponent a is in- 
troduced which allows the transference of the inequality f n (t)> f(t) to 
fn, a (t) ^ f(t). A criterion is given for deciding convexity with exponent 
a in terms of the transform, /(s). 

The degree of precision concept is applied to the approximation se- 
quence in order to obtain a modified sequence, s n (t), which converges 
more rapidly. For sufficiently smooth functions this method is successful. 
The approximation s n (t) consists of a linear combination offo(t), . . . , 
f n (t) or of fo,„(t), .... f n ,„(t) and hence is easily applied. Its efficacy is 
demonstrated in the examples of this paper. Unfortunately the im- 
provement in rapidity of convergence is so strong that the map from f(t) 
to s n (t) is no longer positive, consequently many of the desirable struc- 
tural preservation properties of the L n are lost in favor of greater nu- 
merical accuracy. 

An attempt is made to enhance the rapidity of convergence of \f n (t)} 
while simultaneously retaining the positivity of the map. This is ac- 
complished by the construction of a new sequence, h n (t), which is also 
a linear combination of f (t), . . . , f n (t). As is to be expected, however, 
the improvement is not as great as is realized with the sequence s n (t). 

The pointwise error estimate developed in Part IV may be used as a 
correction device on f n {t) or f n , a (t) to improve further the accuracy of 
computation. This, however, in the absence of an error estimate for the 
modification, must rely on one's understanding of the specific problem 
for ascertaining the reasonableness of the result. 
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An application of the methods of this paper to the renewal function 9 
in the theory of renewal processes is made in Part VI. The remarkable 
accuracy of the simplest of the approximations f (t), fi(t) is notewor- 
thy. 

Part VII presents applications of the techniques to the covariance and 
recovery functions of Erlang blocking models used in telephone traffic 
theory. 11 For the covariance function, the initial approximation, fo, a (t), 
is excellent; however, in the case of the recovery function it was found 
that fo, a (t), h,a(t) might not be considered sufficiently accurate, ac- 
cordingly the linear combination, si(t), was used. This provided suffi- 
cient enhancement of accuracy. 

The generating function, G(z,t), for the f n (t) can sometimes be used 
to obtain an explicit construction of the sequence. Some examples of this 
nature are treated in Part VIII. 

Applications of the methods herein have been made to the comple- 
mentary distributions of waiting time in M/G/l queues. Also B. W. Stuck 
and E. Arthurs have successfully applied these techniques to the study 
of models of computer systems. 

There are questions of an exclusively mathematical character which 
have not been touched upon, e.g., a semigroup interpretation and satu- 
ration phenomena. It is felt that these would be outside the essentially 
practical thrust of the paper. For some theorems which are applicable 
to the operators of this paper see Ref. 5. 

A short table of operations on f{t) and their corresponding maps under 
L n is included to facilitate application of these methods to the con- 
struction of approximations. 

II. WIDDER INVERSION— REPRESENTATIONS 

Let the transform f{s), 

/( s )= f" e - su f{u)du (1) 

exist for s > 0, then 



(-D w 
n\ 

in which 



n+l7(n)( s )=- C e - su U n f(u)du 

n\ Jo 



(2) 



? {n) (s)=£^Hs). (3) 

The function (s n+1 /n\)e- 8U u n is a probability density function on (0, 
oo) for s > 0, n > whose mean is (n + l)/s and variance (n + l)/s 2 . When 
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s = (n + l)/t, the mean and variance are t and t 2 /(n + 1) respectively. 
One has: 1 

Theorem 1 iWidder). Let the transform, /(s), of fit) exist for s > 0, let 
fit) be continuous at t and bounded on [0, »], then 

(-1)" I 

lim — s n+l J^(s) = f(t). 

n-oo Tl\ \s=(n+l)/t 

The convergence is uniform in every finite closed interval throughout 
which f(t) is continuous. 

Proof. Korovkin's theorem on sequences of positive functionals. 2 

The inversion theorem, in the above form, had already been stated 
by Feller 3 who used the law of large numbers to effect the proof. It is the 
purpose of this paper to study the transformation 



L n f = fn=^r-S» + l}(» 



W\ (4) 

n\ \s=(n+\)/t 



so that f n may be effectively used as an approximation to /. The repre- 
sentation of f n directly in terms of / is obtainable from (2); thus, 

fnit)= JJ ' g n it,u)fiu)du (5) 

in + l) n+1 

g n it, U) = e -[(n+l)/t]u u n n > 0. (6) 

Alternative forms which will be found useful are: 

fnit) = 2-J- X" *("' in + ^ f) f{U)dU (?) 

f n it) = in+l) C" f(n,(n+l)u)f(tu)du (8) 

^' o) = e "°?^T) (9) 

in which \pix, a) is the Poisson probability distribution, 

fnit)= r Mn (±)f( U )— (10) 

Jo \u/ u 

M n ix)=^ f e -(n+l)/x x -n-l (n) 

in which the representation is by means of convolution on the half line 
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(Mellin convolution), and 

gM= fS„ KM8iv ~ m U2) 

K n ( v ) = in+ ) )n+1 e-l" +1) *-* +l)e - v (13) 

t = e\u = e*. f(t) = g(v), fn(t) = 8n(v) (14) 

in which the representation is by means of convolution on the whole real 
line (Fourier convolution). 
The conditions of Theorem 1 are relaxed below. 

Theorem 2. Let the transform, J(s), of f(t) exist for s > c, let f(t) be 
continuous at t and let/(t) = 0(e ct )(t — °°); then, 

lim /„(t) =/(*)• 

The convergence is uniform in every finite closed interval throughout 
which f(t) is continuous. 

Proof. The representation (5) may be written as follows: 

f (t) = (" + D" +1 f" e -[(«-m+i)/t]u une -m/t^( u )d w . (15) 
/n n!t n+1 Jo 

For all t in some finite closed interval, m may be chosen so that 

e -(m/t)uf( u ) = 0(1)( W -» oo) 

hence Korovkin's theorem is again applicable and the conclusions fol- 
low. 

III. PROPERTIES OF /„(') 

Jensen's theorem applied to (5) proves 
Theorem 3. f(t) is convex on (0, «) 

^f{t)<f n (t),t>0,n>0. 

The value of an approximation method is greatly enhanced when the 
approximating function preserves the shape of the original and coincides 
closely with its behavior in the neighborhoods of zero and infinity. 
Theorems 4, 5, 6 establish the desired properties. 

Theorem 4. a < f{t) < b =* a < f n (t) < b;a, b arbitrary real. 

Proof. Direct evaluation shows that 

L n f = f,f=<x + 0t. (16) 
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The positivity of L n implies 

f<b=*L n f<L n b = bL n l = b. (17) 

Similarly for the lower bound. 

The derivatives of f n (t) may be related to those of f(t) through use of 
(8); thus 

Theorem 5. Let/< r >(t) be continuous and 0(e ct )(t — °°), then there is 
an m so that/<, r) (£) exists and is continuous for n > m and 

f l n r) (t) = (n + 1) f " rp(n, (n + l)u)u r f(r){tu)du. 
One may set m = if / (r) (£) = 0(l)(t -> »). 

Proo/. For m sufficiently large, the integral of the theorem converges 
uniformly in t; hence the representation (8) may be differentiated under 
the integral sign r times. If/W(t) is bounded then m = is a permissible 
choice since one still has uniform convergence. 

Corollary 1. f<r) > -» /< r) > 0, n > m. 

Proof. This follows from the positivity of the kernel. 

The above corollary implies that if/ has a continuous derivative and 
is monotone then /„ is monotone, and if / has a continuous second de- 
rivative and is convex then f n is convex. A stronger structural result will 
be obtained in Theorem 6. One also has that if/ is completely or abso- 
lutely monotonic then /„ is completely or absolutely monotonic re- 
spectively. 

Corollary 2. /< r) (0+) = A„, r /< r >(0+), n > m, 

T(n + r+ 1) 



*n,r — 



n\(n+ lY 



In particular 



/n(0+)=/(0+) n>m 
/„(0+)=/(0+) n>m. 



Proof. Define A„ >r by 

A n>r =(n+1) f f(n, (n + \)u)u r du 

then evaluation of the integral yields the formula stated. Since the op- 
erator is bounded, the limit statements follow. Also one has X n , = X„,i 
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Corollary 3. Let/ (r) («>) < °°; then 

f { n r) (°>) = X„. r /<'>(») n > 
In particular 

/»(»)=/(») ri>0, 

Proo/. Dominated convergence allows the interchange of limit and in- 
tegral. 

The following concepts will be needed to establish further structural 
properties. 

For an arbitrary sequence in (— °°, °°), — °° < H, t^ . . . , tt < «, the 
number of changes of sign is called the variation of the sequence and will 
be indicated by v(t\, t 2 , . . • , te)\ thus 

o(3, -1,0, 2, -2) = 3 (18) 

u(l,2,4,6) = 0. (19) 

One sets u(0, 0, . . . , 0) = -1. Let/(t) be defined on (0, »), and let < h 
< ti < . . . < t e < °° be an arbitrary, ordered sequence in (0, °°), the 
quantity sup v{f(ti), f(t 2 ), ■■■, f(t e )), in which the supremum is taken 
over all sequences, i.e., for all choices of (ti, t 2 ,---, te) and for all £>l, 
is called the variation of/ and will be indicated by v(f). A transformation 
L on a given class of / will be called variation diminishing if and only if 
u(Lf) <v(f) for every / in its domain. The definition used here is adopted 
from Hirschman and Widder. 4 

Let 0(»7) be a frequency function on (-«, «>), that is, 



0(i/) > 0, f° (j>(v)dr, = 1 



(20) 



and let 



Define E(s) by 



0(8) = £ 



e-"><i>(v)dri. (21) 



E(s) = X*)- 1 . (22) 

Then a theorem of Schoenberg 4 states that the transformation 

Tg = f ~ 4>(S)g(v - 0d£ geBC(-~, 00) (23) 
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is variation diminishing if and only if 

E(s) = e~ Cs2+bs U(l~—) e s/ak (24) 

C>0,b, a k real, £ l/a\ < ». 

The designation gcBC(—°°, °°) means that g(rj) is bounded and contin- 
uous on (— °°, »). It may be observed that the mean of is 6 and the 
variance 

2C + £ l/al 

The Laplace transform of K n (r}) (13), 

K n ( S ) = J^e-^KMdv (25) 



is 



n!(n + IV 
This may be written in the following forms 

*.«- f^ ft (l + f) (27) 

(n + l) s fe=i \ fc/ 

^(s)- 1 = e s " n n (l+l)e~ s/k (28) 

v n - Zrc (n + 1) + t - £ - (29) 

in which 7 = 0.5772157 is Euler's constant. Thus by Schoenberg's the- 
orem, the transformation Tag = g n , g(BC(-<°, °°), defined in (12) is 
variation diminishing. The mean and variance of K n are respectively v n 
as given above, and o\ given by 

2 7T 2 n 1 

^"e""^' (30) 

Since the map t = e v is monotone, the following theorem has been es- 
tablished. 

Theorem 6. The transformation L n defined on ftBC(0, °°) is variation 
diminishing, i.e., 

v(f n )<v(f). 

Corollary. f n does not cross any straight line more often than /. In par- 
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ticular, if/ is monotone then /„ is monotone, and if/ or — / is convex then 
/„ or — /„ is convex. 

Proof. From (16) and Theorem 6, 

u(L n (f -a- fit)) = u(L n f -a-(3t)<v(f-a- (3t). (31) 

Clearly / is monotone <=> v(f — a) < 1 for arbitrary a, and / or — / is convex 
<=» u(f — a — fit) < 2 for arbitrary a, fi. 

It is clear from (8) that the approximating sequence to f(at)(a > 0) 
is f n (at); this may be expressed in a more illuminating way as follows. 
Define the operator A by 

Af(t)=f(at) a>0 (32) 

then A and L n commute; thus 

L n Af = AL n f (33) 

Hence the eigenfunctions of A, which are t r , are also the eigenfunctions 
of L n . In fact one easily obtains 

L n t r = \ nir t r , r>0,n>0. (34) 

It may be observed that if L n is defined by (8) instead of by (4) then (34) 
remains valid even for r < provided n is large enough. 

Other operations with the same eigenfunctions will also commute with 
L n . Of importance in discussing the convergence of fn\t) is the opera- 
tor 

6 = t -f (35) 

at 

One has 

Theorem 7. The operators L n (n > 0), 6 commute; thus 

Ln&f) = 6%. 

It is assumed that the rth derivative of/ exists and is continuous on (0, 
-). 

Proof. It is observed that the eigenfunctions of 6 are t r ; alternatively the 
result follows directly from (8). 

Corollary. Proceeding inductively, one can now establish that if / (r) is 
continuous and 0(e ct )(t —* <») then 

lim/<r> = / (r) 

n— »eo 

In addition to shape preserving properties, another way of assessing 
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the adequacy of an approximation process is by comparison of moments; 
the rth moment of fit) is here taken to be f£t r f(t)dt. The Mellin 
transform, ](s), given by 

J(s)= §~t*-*f(t)dt (36) 

is the appropriate tool. Since the transform of M n (s), eq. (11), is 
^ 0. + OT0.-.+ 1) 

n! 
one has from (10) 

(n+l)T(w-«+l) 7/ ., 

/„(«) = J A*)- (38) 

n! 

For the following, it is convenient to use the factorial symbol 

n (0) - i f n (r) = n ( n - i) . . . ( n - r + 1), r > (integral). (39) 

The following theorem may now be stated. 

Theorem 8. Let the rth moment of f(t) exist (r > 0, integral); then the 
rth moment of f n (t) exists for n > r and one has 

J. a. (n + ~\ ) r+1 /»°° 

Special cases are 

Cf n (t)dt = r}L± ± C"°f{t)dt n>\ 

Jo n Jo 

J(n+ I) 2 /•° 8 
tf n (t)dt= , ' f tf(t)dt n>2. 
o n(n — 1) Jo 

When fot~ l f(t)dt exists, an interesting special case of (38) occurs for 
s = 0; thus 



(40) 



(41) 



C~ t~ l f n (t)dt = fj t-*f(t)dt. 



(42) 



Formulae (40), (41) may be used to ensure equality of moments. Thus 
if it is required that the zeroth order moments agree, then, according to 
(40), one may use as the approximating sequence nf n (t)/(n + 1). If the 
zeroth and first moments are to agree simultaneously, one may use a 
linear combination of /„(*) and f m (t); for example, %fs(0 _ %/2tt). 

Another set of moment relations may be obtained from (7) involving 
sums of f n ((n + l)t). These are given in 
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Theorem 9. Let the rth absolute moment of f(t) exist; then 
£ » w /n((n + DO = t" r_1 f " Wf(u)du. 

Proof. One has 

£ n^Wn, a) = a r (43) 

n = 

and, from (7), 

/„((n + l)t) = t~ l JJ 4, (n, -^j f{u)du. (44) 

Multiplication of both sides of (44) by n (r) and summing, the result 
follows on interchanging summation and integration. Dominated con- 
vergence justifies the interchange. 
Another property of f n (t) as a function of n is given in: 

Theorem 10. Let f(t) > on (0, «), 0(e ct )(t -* °°), then there is an m so 
that 

m + 2\«+ 2 
/ n+1 (£) < e-i (^y) /»(*), t > 0, n > m. 

Proo/. Since 

(n + 2)i^(n + 1, (n + 2)u) = ue~ u (- Y (n + l)\p(n, (n + l)u) 

\n + 1/ 

(45) 
one may write 

fn+i(t) = (n + 1) f " *(n, (n + l)u)ue~ u (^^^Y^fit^du. 
Jo \n + 1/ 

(46) 
Observing that ue~ u < e~ l , the inequality follows. 

Corollary. f{t) > => — — is monotone decreasing in n for all t > 0, n 

n + 1 
> m. 

Proof. One has 

yii^M),.!/^ 1 . (47) 

/i + 2 n + 1 \n + 1/ 
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The result follows since 



n + 2\"+i 



\n+ 1/ 



(48) 



A stronger monotonicity property of f n (t) is stated in Theorem 21. 

A useful recurrence relation allowing one to compute the members 
of the sequence f n (t) successively starting with f Q (t) is given in the fol- 
lowing theorem. 

Theorem 11. 

Proof. Define f n (s) by 

/n(s) = -^T^ n+1 / (n) (s) (49) 

n! 



then, by (4), 

fn(t)=f n (s) 



s-(n+l)/( 



(50) 



One has 



4/.W - (n + 1)/»W - (n + l)/ n +i(s) (51) 

as 

fn+i(«) -/»(•) — r^T-/" 00 - (52) 

n + 1 as 



Thus 



fn+l(t) = \fn(s) ~ 4t S T^ (S) 1 

L n+l ds Js=(n+5 



(53) 

•2)/t 

The recurrence relation is now obtained on performing the substitution 
for s. 

A useful alternative method of presenting the structure of the entire 
sequence \f n \o in terms of /o is by means of a generating function. This 
is given in the theorem below. 

Theorem 12. f(t) is bounded on (0, ») => 

£ z»f„Un+ 1)0 =7^-/0 (t-M- 

n=0 l -2 M-2/ 

The series is convergent for \z | < 1 and analytically continuable in the 
half plane Re z < 1. 
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Proof. The formula follows from (7) after interchange of summation and 
integration. The series is clearly convergent for \z | < 1 while dominated 
convergence justifies the interchange for Re z < 1. 

The case r = of Theorem 9 provides the following corollary. 

Corollary. 

f(t)eL(0, »)=» lim -i-/ (-M = f-l f"/(a)da. 

2-.1- 1 — 2 \1 — 2/ »/0 

The Mellin transform, /(s), of f{t) may be directly obtained from its 
Laplace transform, /(s), by use, for example, of (38) for n = 0; thus 

7/ \ /OOP ,r.. x 

m= ^T7V (54) 

Accordingly one may now write (38) in the form 

/n(«) = T=t; /o(s) (55) 

n!r(l -s) 

or, equivalently, 

] n (s) = (n+l)°( n ~ S yfo(s). (56) 

At times (55) or (56) provides a convenient alternative to Theorems 11 
and 12 when f n (t) is required as a function of n. 

The range of applicability of the Jensen inequality of Theorem 3 may 
be extended by use of Mellin or Fourier convolution. A sequence \f n (0}«=o 
will be called an approximation sequence if there is an f{t) so that/„ (t) 
= L n f(t). Let * designate Mellin convolution; then 

Theorem 13. f n *g is the approximation sequence for f*g. 
Proof. One has 



L n (f*g) = M n fg = L n fg (57) 

thus, 

L n (f*g) = (L n f)*g = f n *g. (58) 

The converse of Theorem 3 is also true. 

Theorem 14. f n (t) > f(t) for all n > 0, t > 0, f(t) is bounded on (0, «) -> 
f(t) is convex on (0, °°). 

Proof. The result follows from Theorem 8 of Karlin and Ziegler. 5 

Corollary, f convex on (0, °°), g > on (0, °°) ==> f*g convex on (0, »). 

682 THE BELL SYSTEM TECHNICAL JOURNAL, MARCH 1978 



Proof. One has from Theorem 3 

fn>f (59) 

and, since g > 0, 

fn*g^f*g- (60) 

Since, by Theorem 13, f n *g is the approximation sequence to f*g, ap- 
plication of Theorem 14 proves the corollary. 

It may be observed that the inequality of (60) remains valid when * 
is interpreted as Fourier convolution although, in this case, f n *g is not 
the approximation sequence for f*g. 

Another set of convexity results may be obtained from (8) by consid- 
ering logarithmic convexity. 

Theorem 15. If f(t) is log-convex on t > 0, then f n {t) is log-convex for 
n > 0, t > 0. 

Proof. Equation (8) and the additivity of log-convex functions. 6 
Further one may state the following inequalities. 

Theorem 16. If f(t) is log-convex on t > then 

f{t)<e L ^m<f n { t ). 

Proof. The inequality on the left follows from Theorem 3 applied to 
£nf(t); the one on the right is a consequence of the geometric mean- 
arithmetic mean inequality. 

IV. ERROR ESTIMATION 

Error estimates take different forms depending on the class of func- 
tions for which they are intended and whether or not they are bounds 
or pointwise estimates. From a practical point of view the pointwise es- 
timate is the most useful provided it may be easily evaluated in terms 
of the approximation itself. The next three theorems provide error 
bounds for different function classes; the fourth theorem provides an 
approximate formula for the pointwise evaluation of error, while (112) 
does the same but in terms of f n (t). The error of approximation, t n (t;f), 
is defined by 

<n(t;f) = fn(t) ~ fit) (61) 

Theorem 17. Let f(t) be continuous on (0, °°); then 

Mt;/)|^-7=sup|/(t)|. 
vn + 1 t>o 
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Proof. One has 

<n(t;f) = fjsnit, u)\f(u) - f(t)\du (62) 

\e n (t;f)\ < y g n (t,u)\f(u)-f(t)\du (63) 

h„(i;/)| < sup |/(i)|- f" g n (t,u)\u-t\du (64) 

t>o Jo 

Un(t;f)\ < sup 1/(01 . ( C~g n {t, u)(u - t) 2 du U (65) 

t>o I Jo 

\e n (t;f)\ <— ~sup\f(t)\. (66) 

Vrc + 1 t>0 

The last inequality follows because the mean and variance of g n (t, u) 
are t and t 2 /(n + 1) respectively. 

Theorem 18. Let f(t) be continuous on (0, °°); then 

|6„a;/)|<-^-sup|/(0|. 
In + Z t>o 

Proof. The Taylor expansion of f(u) about t has the form 

f{u) = f(t) + (t- u)f(t) + \(t- u) 2 'f(0 (67) 

in which £ lies between £ and u. Thus 

e " U;/) = 2^+T fe), ^ (0,oo) - (68) 

The inequality of the theorem now follows. 

The next theorem provides an error bound which is uniform for te[0, 
»]. For this purpose the absolute first moment of K n (ri) (13), a n , is 
needed; thus 



a n = J_ a> K n(v)\v\dri. 

Theorem 19. Let f{t) be continuous on (0, °°); then 



(69) 



| «„(*;/) | <a n sup \tf(t)\. 

t>o 



Proof. One has from (12) 

e„(e";/) = J " K n ( v - S)\g(S) - g( v )\d$ 
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(70) 



|e„(e»;/)|< sup \g'( V )\ • f" K n ( v - 0\v - t\dl- (71) 
| € „(t;/)| <a„sup|t/(0|. 

t>0 

Corollary. The convergence of f n (t) to /(£) (n -» °°) is uniform for £e[0, 

Proof. It is necessary to show that 

lim a n = 0. 

n— »oo 

The expression for «„ (69) is rewritten as follows 

an = Pn f_l K ^) n+1 \v\dv (72) 

1 in + 1\«+1 

p n = -r( ) .IfW-e 1 "^. (73) 

n\\ e I 

Use of the power series expansion for e~ v yields 

«n~Pn r e-^+W\ v \d v = — (? L±1 ) n . (74) 

J-~ en! V e / 

Stirling's formula now shows that 



7rn 



(75) 



Some numerical values of a„ are «o = 1.0160, a\ = 0.6388, c*2 = 0.5006, 
a 3 = 0.4247, a 4 = 0.3751, a 5 = 0.3396, a 6 = 0.3126, « 7 = 0.2911; the as- 
ymptotic formula (75) is sufficiently accurate for n > 7. 

To continue the study of e n (£;/), it is useful to obtain an explicit for- 
mula of Peano type, that is an integral transform of/. 
Let 

x+ = x x>0 (76) 

= x < (77) 

then the Taylor expansion of f(t) with remainder is 

fit) = /(0) + W)t + C " (t - v) + f{v)dv. (78) 

From (5) and (16), one has 

fn(t) = f(0) + f(0)t + C~Hv)du C~g n (t,uKu-v)du. (79) 

JO Ju 
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Thus 

*»(«,/)- fj E n (t, v)f(v)dv (80) 

E n (t, V) = fjgnit, V)(U - V)du - (t - v)+. (81) 

The kernel E n (t, u) (the Peano kernel for error representation) is, 
clearly, 

E n (t, u) = L n (t - u) + -(t- u)+. (82) 

The explicit evaluation of the kernel may be most simply carried out by 
means of (24) since the Laplace transform of (t — u)+ is e~ sv /s 2 . Let 

Sn(x) =±^ (83) 

7 = jl 

and 

£ n (a) = e-l" + » a [S n «n + l)a) - aSn-din + l)a)] (84) 

then 

L n (t - v)+ = t£ n (a) a = v/t (85) 

E n (t,v) = t[£ n (a)-(l-a) + ]. (86) 

In particular, one has 

E (t,v) = te-» /t -{t-v) + (87) 

E&, v) = (t + v)e~ 2v ' t -(t- v)+. (88) 

Since (t — u)+ is a convex function of t for each u, (82) and Theorem 
3 establish 

E n (t, v)>0 for all t > 0, v > 0. (89) 

The moments of the kernel, E n (t, v), may be obtained by substituting 
the functions f(t) = t r (r > 2) into (80), and using (34) and (61) for 
evaluation of e n (t;t r ); the following is obtained: 



r+2 r > 0. (90) 



Jo (r+l)(r+2) 

In particular 

C~E n (t,v)dv = t -^^— (91) 

Jo 2 n + 1 

rvE n (t,v)dv=^^^-- (92) 

«/o 6 (n + 1)^ 

One may now obtain an approximate evaluation of e n (t;f). 
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Theorem 20. Let f(t) be continuous on (0, «) and 0(e ct )(t — °°); then 
there is an m so that 

, - t 2 .. / 3rc + 5\ 

also if /(£) is convex, then the approximation is a lower bound. 

Proof. The one point Gaussian quadrature formula for foE n (t, u)f(u)du 
is of the form Af(a) in which the constants, A, a, are determined by re- 
quiring the quadrature to be exact for all linear functions. Use of (91), 
(92) now yields the formula of the theorem. The inequality follows from 
the nonnegativity ofE n (t, u) (89) and Jensen's inequality. 

Since by the Corollary to Theorem 7, f n (t) approximates f{t), in 
practice the required value of f(t) is approximated either from the an- 
alytic form of f n (t) or numerically from a table or curve already com- 
puted for /„(£)• 

At this point another property of the sequence {/ n (t)}o can be 
proved. 

Theorem 21. Let f(t) > 0, continuous on (0, °°), and 0(e ct )(t -* «>); then 
there is an m so that 

fn+i(t) < fn(t) for all t > 0, n > m. 

Proof. Clearly the monotonic decreasing character of /„ (£) as a function 
of n will hold if e n (t;f) has the same property. The nonnegativity of f(t ) 
and (80) shows that the result is implied if E n (t, u) is monotonically 
decreasing in n; in turn, by (86), this will follow if £ n (a) is monotonically 
decreasing in n for each a > 0. From (84), by direct calculation, 

-f £ n (a) = -e-^ + V°S n ((n + l)a) (93) 

da 

d 2 (n + l) n+1 

-J~Jn(a) = { / a"e-<"-"K (94) 

da z n\ 



Let 



then, from (94), 



h n {a) = e n - l {a)-£ n {a) n>\ (95) 



d 2 d 2 

-r^h n (a) = r n (a)—-£ n - l (a) (96) 

da 1 da 2 

I l\ n+1 
r„(a) = l-(l + -J ae~°. (97) 



INVERSION TECHNIQUE FOR LAPLACE TRANSFORM 687 



It is clear from (94) that the sign of 

— Ma) 
da 1 

is the same as that of r n (a). There exist two points < a (n) <ai(n) with 
the following properties: 

r n (a) > < a < a (n) (98) 

r n (a) < a {n) <a< cu(n) (99) 

r n (a) > a > <n(n). (100) 



Since 



it follows that 



*„(0) = 1,-7-^(0) = -1, n>0 (101) 

da 



h n (0) = 0,4-h n (0) = 0, n>l. (102) 

da 

One has the following integral representations for h n (a): 

hnia) = So db X r " (c) ^i^- i(c)dc ' (io3) 

/»»(a)= f"d6 f r B (c)j74-i(c)rfc. (104) 

Ja Jb dC A 

Thus (98) and (103) imply 

h n (a)>0 0<a<a (n); (105) 

similarly (100) and (104) imply 

h n (a) > a > ai(re). (106) 

The function h n {a) cannot be negative in (ao(n), a\{n)) since then it 

would have at least one local minimum; however, (99) shows that in 

(a (n),ai(n)) 

d 2 
— h n (a)<0 

which is a contradiction. Thus 

h n (a)>0,a >0,n>l (107) 

and the theorem is proved. 

It is possible to estimate conveniently e n {t;f) directly from/„(£) if f n (t) 
and /„ (t) are readily obtainable, at least possibly numerically from values 
offnit). From (38) one has 
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f(s) = M n (s)-%(s). (108) 

Expansion of T(n + 1 — s) into a power series in s and substitution into 
(37) provides the following series 

M„(8) = l + V n S + ^"^S 2 + ... (109) 

M n (s)- 1 = 1 - » n s + " n ~ an s 2 + ... (110) 

Thus 

<n(s;f) = [p n s + £ ^f J ^s 2 + • . ] In (HI) 

2 2 

*„(*;/) « -r n Of n {t)+^—^0*f n (t). (112) 

To facilitate the use of (112) some values of the coefficients are given 
in Table I. 

The following readily obtained asymptotic formulas may be used for 
values of n beyond the table: 

+ 7T. - ... (113) 



2(n + 1) I2(n + l) 2 

1 + x 



" n + 1 2(n + l) 5 



+ 



2 2(n + 1) 8(n + l) 2 

V. ENHANCEMENT OF ACCURACY 

The excellent behavior of the operator L n in constructing approxi- 
mations to a given f(t) which preserve its structural properties and its 
limiting values and which provide inequalities exacts a penalty in the 
form of slow convergence. A high value of n is required to attain high 

Table I — Coefficients 



n 


»n 


«l 


(<* " "2)/2 





0.5772 


1.6449 


0.6559 


1 


0.2704 


0.6449 


0.2859 


2 


0.1758 


0.3949 


0.1820 


3 


0.1302 


0.2838 


0.1334 


4 


0.1033 


0.2213 


0.1053 


5 


0.0856 


0.1813 


0.0870 


6 


0.0731 


0.1535 


0.0741 


7 


0.0638 


0.1331 


0.0645 


8 


0.0566 


0.1175 


0.0572 
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numerical accuracy. In many practical problems, fortunately, very high 
accuracy is not needed; notwithstanding, the value of n required may 
still be inconveniently high. Considering that one starts with the Laplace 
transform, J(s), of f(t) and uses (4), or constructs f (t) and uses the re- 
cursion of Theorem 11, a high value of n implies obtaining a corre- 
spondingly high order of derivative of J(s) or of fo(t) which can be a 
time-consuming operation. Thus it would be useful to modify the basic 
approximation, L n f, while still preserving many of its original charac- 
teristics so that the accuracy for a given value of n may be increased. 
In many cases the transform, /(s), has the property that for some a 
> 0, f(s - a) converges for s > 0. This property is used to construct a new 
approximation, f n , a (t), defined by 

UAt) = e-«*L n (e«W)) (H4) 

and, correspondingly, a new operator 

L n ,af = fn, a . (116) 

The following theorem permits the construction of /„,„(*) directly from 
fnit). 

Theorem 22. 



fn.M = fn . 

/ at \n+i I at 

\ n + 1/ v " n + 1 

Proof. From (5) and (6), one has 
(n+l) n+1 



Lnf- t L f" e-^ n+1 ^ u u n f(u)du (116) 

n\t n+1 Jo 

r (pa tf) = (n+ 1)n+1 f~ e -[(.n+l)/t)u+au u nf( u )du. (117) 

nK n n\t n+1 Jo 



Thus 



L n f\ l- a t/(n + l) 

= /l-^^-V +1 (n + 1)fl+1 r" e -[("+i)Alu+«« M »/(a)du. (118) 
V n + 1/ n\t n+l Jo 

Comparison of (118) with (117) shows that 

L n (e°'f) = /„ I I (119) 

\ n + 1/ x n + 1 

hence, the result follows from (114). 
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The approximations f n , a (t) satisfy theorems similar to those proved 
for f n (t); however, modifications are required. Only Theorem 3 will be 
discussed. A function f(t) will be said to be convex over an interval with 
exponent a if and only if e at f(t) is convex over the same interval. 

Theorem 23. If f(t) is convex on (0, <») with exponent a then 

f(t)<f n>a (t). 

Proof. One has from Theorem 3, 

e«<f(t) < L„(e at f(t)). (120) 

The result now follows from the nonvanishing of e at and (114). 

The error of approximation by /„,„(£) will be designated e n .a(t;f) and 
defined by 

tn.aitj) = fn.a(t) ~ fit). (121) 

Clearly 

*n,«(t;f) - e- at e n (t;e at f). (122) 

also, if the condition of Theorem 21 is satisfied, one has 

€„,«(*;/) >0. (123) 

One of the useful aspects of the approximation, /„,„ (t), is that it more 
accurately reflects the asymptotic behavior of f(t)(t -* «> ) than f n {t) does 
for a given value of n. In the later applications this will be an important 
characteristic. 

Clearly, ordinary convexity corresponds to convexity with exponent 
zero; however, the following theorem relates convexity with exponent 
a to log-convexity. 

Theorem 24. Let f{t) be continuous on some interval I; then f{t) is log- 
convex on / if and only if it is convex with exponent a on I for all a. 

Proof. One has 

£n (e" l f(t)) = at + £nf(t) (124) 

hence e at f(t) is convex with exponent a on J for all a if f(t) is log-convex 
on /. The derivative condition for convexity with exponent a on / is 

f(t) + 2af(t) + a 2 f(t) > on / (125) 

and the derivative condition for log-convexity on / is 

f(t)f(t) - f(t) 2 >0 on I. (126) 

The choice 
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which is always possible since f(t) > on J, in (125) verifies (126). 

Convexity with exponent a and, hence, by Theorem 24, log-convexity 
may be decided by means of the Laplace transform and the use of the 
Hausdorff-Bernstein theorem. 7 

Theorem 25. Let f(t) be continuous on (0, °°), then f(t) is convex with 
exponent a on (0, °°) if and only if 

(s + a) 2 J(s) - (s + 2a)/(0+) - /(0+) 

is completely monotonic in s on (0, «>) and is absolutely convergent on 
s>0. 

Proof. The expression cited is the Laplace transform of 

e - a tJlL {e *t f{t)) 
at' 

whose nonnegativity is the necessary and sufficient condition for con- 
vexity of f(t) with exponent a. The Hausdorff-Bernstein theorem now 
completes the proof. 

It may be observed that the quantities /(0+), /(0+) are obtainable 
from 

lims/(s)=/(0+) (128) 

lim \s 2 f(s) - s/(0+)} = /(0+). (129) 

S -»OD 

Another method of enhancement is related to the concept of "degree 
of precision." An approximation operator T, i.e., Tf « /, in which the 
functions t r , suitably restricted to an appropriate interval (r > 0, inte- 
gral), are in its domain, is said to have degree of precision k if Tf = t r 
for < r < k and Tt r ^ t r forr = k + 1. Thus the singular operators L n 
studied here have degree of precision one. 

The enhancement method consists of the following: coefficients 5/(0 
< j < k — 1) are determined by the moment conditions 

Z6jLj(t r ) = t r 0<r<k (130) 

;'-o 

and, accordingly, the linear combination 

s k -i(t) = k £f>jfj(t) (131) 

j-o 

is now taken to approximate f(t). Clearly the map from / to Sk-i has 
degree of precision k, however, unfortunately, it is not positive. If / is 
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sufficiently smooth there will result a significant improvement in ac- 
curacy over the use of fk-\ alone. The system (130) may be expressed 
in terms of \ n>r as follows 

k-i 

£ 8j\ jr = 1 l<r <k. (132) 

/-o 

Accordingly special cases of (131) are 

si = -/o + 2/x (133) 

S2 = ^o-4/i + ^/ 2 . (134) 

A method of enhancement consisting of a linear combination of the 
f n (t) similar to (131) which, however, retains the positivity of the map 
will now be constructed. The accuracy attained will usually not be as 
great as that of (131) for a given set of values \fj(t))%. The new sequence 
will be designated h n (t) and is defined by 



h n (t) = E Pjfj(t). (135) 

1=0 



Define W n (u) by 



W n (u) = E Pjti + WO', 0" + Dm) (136) 



then the coefficients, pj, are constrained by 



E P j = 1, W n (u) > for all u > 0. (137) 

7 = 



Theorem 26. 

11 



\h n (t)-f(t)\ <t sup |/(t)| 

t>0 



- ■ 1/2 

;?o7+ i 



Proof. From (8), one has 

h n (t) = JJ W n (u)f(tu)du. (138) 

Also, from (137), 

f" W/ n (w)du= 1 



(139) 



hence 



h n {t)-f(t) = J^° W/ n (u)[/(tu) -/(t)]d«i. (140) 
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The nonnegativity of W n {u) now permits the following inequality 

\h n (t)-f(t)\ < r°W n (u)\f(tu)-f(t)\du; (141) 

hence, 

\h n (t)-f(t)\<t sup \f(t)\ C"w n (u)\u-l\du. (142) 

t>o Jo 

The Cauchy-Schwartz inequality applied to (142) yields 

\h n (t)-f(t)\ <t sup |/(t)| ( f " W n (u){u - l) 2 du 
t>o (Jo 



11/2 

• (143) 



Evaluation of the integral in (143) provides the inequality of the theo- 
rem. 

Consider the sum, S, defined by 

S = L t^- (144) 

j=oj + 1 

then, in order to obtain the best approximation, the pj must be chosen 
to minimize S besides satisfying the conditions of (137). One has 

n (j + IV'+l . . 

e»W n (u) = E ., pje-'W. (145) 

j=o Jl 

Let 

2 = e l ~ u u (146) 

then the constraint of (137) may be written 

n (j+ iy+i . . 

E - — 7, e-J Pj zJ > < z < 1. (147) 

j=o Jl 

Define the polynomials P(x) by (147) with z = (x + l)/2; then one 
has 

PW=|/[|^^(*)<2e)-^] (148, 

and 

P(x) > -1 < x < 1. (149) 

The cosine polynomial P(cos 6) is now obtained and written in the 
Fourier form 

1 » 

P(cos 0) = - a + E a ; cos ;'fl (150) 

2 ;=1 
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-it 



in which the ay are obtained from 

P(cos0)cos;'0d0. (151) 

The nonnegativity of P(cos 0) implies the following representation 8 

P(cos0) = |M0)| 2 (152) 

h(d) = x + x 1 e io + ... + x n e in0 (153) 

in which the coefficients x , . . . , x n are real; thus 

aj = 2 "-£ x t x v+j 0<j<n. (154) 

k=0 

When (154) is solved for the pj in terms of x , . . . , x n , the problem of 
minimizing 5 subject to p + • • ■ + p n = 1 becomes that of minimizing 
a quadratic form relative to another quadratic form. 
The optimum pj have been obtained for the case n = 2; the result is 

h 2 {t) = 0.146993/oU) - 0.944260/i(*) + 1.797267/ 2 (0 (155) 
with S = 0.273952. Thus, from Theorem 26, 

\h 2 (t) - f(t)\ < 0.523404i sup |/(t)|. (156) 

t>o 

The estimate of t n (t\f) in (112) may be effectively used to reduce error. 
One may take as an approximation to f(t) the following 

2 2 

fit) m f n (t) + v n 9f n {t) - ^—^ 0%(t). (157) 

In order to improve f n , a (t), the approximate calculation of c n , a (t;f) 
proceeds by use of (122). The practical use of (112), (157) uses difference 
quotients to evaluate Of n (t), d 2 f n (t) from the values already obtained 
for f n (t). Thus, let h > be the distance between consecutive values of 
t for which f n (t) is calculated; then 

y.w -« '-<*+*> -'■<*-*> (158) 

f/.(0«W + ^-" + fc) -W + /-"-*» . (159) 

h l 

The following comment should prove useful in reduction of error. If 
a function g(t) is known which approximates f(t), for example, the 
leading term of an asymptotic expansion for f(t), then one may use 

f(t) m g( t ) + L n (f - g). (160) 
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Evidently an appropriate g(t) should always be sought before con- 
structing practical approximations to f(t). 

VI. THE RENEWAL FUNCTION 

In this section some of the preceding theory will be applied to ob- 
taining approximations for the renewal function, M(t), of a renewal 
stream. 9 Let A(t), with A(0+) = 0, be the interarrival time distribution 
and A(s), given by 

A{s)= C~ e- st dA(t) (161) 

its Laplace-Stieltjes transform, then 

ww-Vttt (162) 

s 1 - A(s) 

The sequence of approximations, M n (t), may now be constructed 
from 



M " (t) = T^klT)- 1 - (163) 



In particular one has 



.,,., 1 , 2 A'(2/t) ,_., 

M * t)m T=m)- l "ii-M2/w (164) 



in which 



A'(s)=-f A(s). (165) 

as 

Let A be the arrival rate, and a 2 the variance of interarrival time, that 
is, 

X-i= C~tdA(t) (166) 

a 2 = C~ t 2 dA(t)-\~ 2 (167) 

then evaluation of the contribution of M(s) at s = provides the 
term 

Thus one may introduce a new function, f(t), by 

M(t) = Xt + a2X2 ~ l + f(t) (169) 

La 
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with 



/0V l-A(l/t) 2 



Since linear functions are invariants of the operators L n , there is no re- 
duction of error when approximating f(t) by f n (t) over approximating 
M (t) by M n (t); however, often f(t) is exponentially dominated and the 
enhancement technique of Theorem 22 is applicable. 
The following example will be considered: 

A(t) = erf \A Ms)= /T 1— (171) 

2 v 1 + 2s 

^ (s) = i ^T^^~T' (172) 

s v 1 + 2s — 1 
Thus, 

M ° ( ' ) = VTT^-x U73) 

MM = VT+17F - i + 7 VT+47t< vTT^-Ip (174) 

Since X = 1, a 2 = 2, one has 

1 2 1 1 

flit) ~ Vl + Alt -1 + t Vl + 4A(Vl + 4/t - l) 2 ~ ^ ~ 2 " 

(175) 

The a transformation of Theorem 22 may be applied to fi(t). As- 
suming fi(t) to be ultimately of one sign, the singularity farthest to the 
right of /(s), namely — V2, coincides with the abscissa of convergence; 
hence, a = %, Table II compares the approximations for M(t) given by 
M\(t), the enhancement procedure of (133), and t + V2 + /1, 1/2(0 with 
more accurate values obtained from the exact solution 

M(t) = f; C e- u u n - l ' 2 du. (176) 

n=0_ / 1\ JO 

V 



r (n + ' ) 



Since M(t) » t + %(t -*• »), the accuracy increases with increasing 
£. This is characteristic of the applications to the renewal function. 

An example will now be considered in which the interarrival time 
distribution is a mixture of exponentials since this is of frequent practical 
use. Accordingly let 
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Table II — Comparison of approximations 



Mi(t) 



Slit) 



t+f+/li(*) 



M(t) 


















0.5 


0.83333 


0.85764 


0.84297 


0.86007 


1 


1.39443 


1.42283 


1.41014 


1.42466 


2 


2.44338 


2.47255 


2.46303 


2.47161 


5 


5.48142 


5.50480 


5.49568 


5.49718 


10 


10.49350 


10.50977 


10.49980 


10.49989 



A(t)= l- — e- t -—e- 2t . 
10 10 



Then, 



A(s) = 



1 + 6 



1 



M(s) = 



10s + 1 10s + 2 

1 20 + 13s 



s 2 17 + 10s 



Also one has 



20 91 

17 289 ' 



fo(t) = - 



210 



289 17t + 10 
Application of Theorem 12 provides 

8400 1 



hit) = - 



flit) = - 



289 (17* + 20) 2 
5.67 X 10 5 1 



289 {lit + 30) 3 
The exact solution for this simple example is 
20 21 



M(t) = — t + — (l - e -U7/io)A . 
17 289 V / 



(177) 

(178) 
(179) 

(180) 
(181) 

(182) 
(183) 

(184) 



Table III compares calculations from M 2(0 and the enhancement pro- 
cedures of (134) and (155) with the exact value. The a enhancement 
procedure with a = 1.7 was not used because it produces the exact re- 
sult. 

This example shows the operation of the enhancement procedures (134) 
and (155); clearly, S2(t) is very accurate since the constraint of positivity 
of the approximation operators is discarded in its construction. 
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Table III — Comparison of calculations 



M 2 (t) 



s 2 (t) 



h 2 (t) 



M(t). 


















0.5 


0.62652 


0.62967 


0.62712 


0.62984 


1 


1.23024 


1.23559 


1.23127 


1.23586 


2 


2.41812 


2.42353 


2.41913 


2.42318 


5 


5.95373 


5.95595 


5.95407 


5.95500 


10 


11.83713 


11.83747 


11.83710 


11.83737 



VII. THE COVARIANCE AND RECOVERY FUNCTIONS 

The study of errors in switch count and continuous scan observational 
methods in telephone traffic engineering is facilitated by use of the co- 
variance function of the number of busy trunks in the Erlang blocking 
system. 11 Specifically let x t be the number of trunks busy at time t in 
an equilibrium MIMIC blocking system with unit mean holding time 
and offered load of a erlangs, then the covariance function, R(t), is 



R(t) = E(xox t ) - (Ex ) 2 . 



(185) 



In order to express the Laplace transform, R(s), it is necessary to in- 
troduce the Poisson-Charlier polynomials 10 - 13 which may be obtained 
from 



Gj(x,a)= ± (-lV-'( ; )*/!a -*(*). 
„=o \v/ \v/ 

They satisfy the following recurrence 

G j+l (x, a) = X ~ J ~ a Gj(x, a) -^Gy-ifo a) 



(186) 



(187) 



G (x,a) = l 



G^x.a) = -- 1. 
a 



Also needed is the function «,(*, a) given by 

Gj-i(x, a) 

otiix, a) = —t 

1 Gj(x,a) 

which satisfies the first order recurrence 



otj+i(x, a) ! = a;(x,a) 

a a 



a l (x,a) = (--l) 



-l 



(188) 



(189) 



The zeros of Gj(x, a) are all positive and simple; in particular, the zeros 
of Gj(—s — 1, a) as a function of s are designated r, and ordered by 
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r j <r j - l <... <n<0. (190) 

In the approximation to be developed, r x will be the dominant root. 
The Erlang loss function, 10 B(c,a), given by 

a c i c aJ 
c\l j=oj\ 

gives the probability that all servers are busy. In the formulae below it 
will be designated simply by B. The mean number of servers busy, n, 
is 

fi = a{l-B) (191) 

and the variance, a 2 , of the number of busy servers is 

<j 2 = n-a{c-ti)B. (192) 

The Laplace transform, R{s), and the covariance, R(t), are 11 
p/ \_ a2 + V 2 , a » acB 

KKS> 1+S 8(1+8) (1+s) 2 



s s(l + «) z 



R(t)= t V o< ^..^-^(L-y 



(194) 



The approximation Ro(t) is 
* o(0 " 1 + t + l + t M 



acfii acfit 2 



+ hS-(-7" 1 -)- < 195 » 



(1 + t) 2 (1 + i) 2 

Since the dominant root is r\ one may choose a satisfying < a < — r\ 
to obtain 

Ro, a (t) = a*e~ at g{t) 

s(t)=-—— - + 



l + (l-«)£ (1 - at)(l + (1 - a)t) 

H 2 /a 2 acBt/a 2 



I -at (l + (l-a)t) 2 

acBt 2 /a 2 / 1 \ / 1QKU x 

H -a c ( — -+ a — 1, a )• (195b) 

(1 - at)(l + (1 - a)t) 2 \ t I 

It is known that the zeros r, are separated by at least one so that 1 — l/(r/ 
- n) > 0, and hence Aj is positive for each ;'; thus R(t) is log-convex. 
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Accordingly, the following inequality is valid (Theorem 23): 

R(t)<R 0ia (t). (196) 

In order to facilitate the use of Ro,„(t) an accurate upper bound for 
ri is needed to provide a suitable choice for a. Such a bound is available 
in Ref. 12. Thus let 

c 1 
£ = £ -c^a~" (197) 

m = £ 2 -2 t-c M a-""J:- (198) 

<-=2 V j=ij 

then 

ri< ° - -1. (199) 

£ + V(c-l)(cm-£ 2 ) 

To illustrate the practical performance of (195b) and (199), calcula- 
tions were made for the cases a = 4, 8, 12 and c = 8 corresponding to 
medium, heavy, and very heavy loads respectively. The corresponding 
equilibrium blocking probabilities are B(8, 4) = 0.030420, 5(8, 8) = 
0.235570, B(8, 12) = 0.422655. Table IV compares the exact and ap- 
proximate values. Figures 1(a), 1(b), and 1(c) compare the corresponding 
curves. 

Table IV — Comparison of exact and approximate values 







a =4 




a = 8 


a 


= 12 


t 


R(t) 


Ro.a(t) 


R(t) 


Ro.a(t) 


T{t) 


Ro.aU) 





3.377 


3.377 


2.564 


2.564 


1.492 


1.492 


0.4 


2.143 


2.145 


1.075 


1.091 


0.312 


0.331 


0.8 


1.365 


1.367 


0.474 


0.483 


0.075 


0.079 


1.2 


0.870 


0.872 


0.212 


0.216 


0.019 


0.020 


1.6 


0.555 


0.556 


0.095 


0.097 


0.005 


0.005 


2.0 


0.354 


0.355 


0.043 


0.043 


0.001 


0.001 


2.4 


0.226 


0.227 


0.019 


0.019 






2.8 


0.144 


0.145 


0.009 


0.009 






3.2 


0.092 


0.092 


0.004 


0.004 






3.6 


0.059 


0.059 


0.002 


0.002 






4.0 


0.038 


0.038 


0.001 


0.001 







The quality of approximation of (199) may be seen from the following 
values of a used in (195b) compared to the exact r x values. 

a — n a 

4 1.1218 1.1215 

8 2.0000 1.9730 

12 3.4778 3.3415 

The transition probabilities P»/(t) — the probability; trunks are busy 
at time t given i trunks are busy at time zero — may all be obtained from 
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the transition probability P cc (t) n ; this probability as a function of time 
is called the recovery function. It may be used in a similar manner to the 
co variance function, R(t), for the study of errors in scan measurement 
techniques 14 ; additionally it is especially important in the analysis of 
telephone retrial models. 

The Laplace transform, P cc (s), and the recovery function, P cc (t), 



are 



Pcc(s) = " + — «c(-s - 1. a) (200) 

s as 

P cc (t) = B- ZBje'Jt (201) 



1 rj ijtj \ n - rj 



(202) 



As for the covariance function, B = B(c, a), and ry(l < j < c) are the roots 
of G c (— s - 1, a) as a function of s. 

In order to apply the a enhancement procedure, the function 

f(t)=P cc (t)-B (203) 

is considered whose Laplace transform is 



f(s) = -\ 1 - B + - a c (-s - l,a) 1 
s \_ a J 



(204) 



It may be observed from (201) and (203) that fit) is log-convex, hence 
the approximations obtained will constitute upper bounds. In order to 
demonstrate the operation of the approximations, the functions /<>,«(* )• 
f ha {t), and 

Sl (t) = 2f ha (t) - fo.At) (205) 

were constructed; they are 

^)-S[ 1 - B+ S*("7 + "- 1 - 8 )]' (206) 
'^)-T^[ 1 -* + i-(-7 + «- 1 -*)] 

\ 2/ 



e- at 2c , / 2 \ 

a tat \ t / 



(207) 



X ~T 



The prime on a c (x, a) indicates differentiation with respect to x. The 
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following recurrence relation is obtained from (189); 

a'j+i(x, a) = - \ja'j(x, a) - l]a;+i(x, a) 2 

a\{x, a) = - - ai(x, a) 2 . 



(208) 



Since si(t) does not correspond to a positive operator, it does not provide 
a bound for P cc (t); one has, however, 

P cc (t) <B + f ha (t) <B + fo, a (t). (209) 

The first inequality follows from Theorem 23 and the second inequality 
from Theorem 21. 

The same cases as for the covariance function were treated. Tables 
V, VI, and VII compare the exact and approximate values. Figures 2(a), 
2(b), and 2(c) compare the corresponding curves. 



Table V — a = 4 



t 


Pcc(t) 


B + f ,a(t) 


B + hmit) 


si(t) 





1.0000 


1.0000 


1.0000 


1.0000 


0.1 


0.5178 


0.5907 


0.5597 


0.5287 


0.2 


0.3304 


0.4280 


0.3856 


0.3432 


0.3 


0.2380 


0.3335 


0.2901 


0.2468 


0.4 


0.1844 


0.2703 


0.2296 


0.1889 


0.5 


0.1497 


0.2249 


0.1880 


0.1511 


0.6 


0.1256 


0.1907 


0.1578 


0.1248 


0.7 


0.1080 


0.1641 


0.1350 


0.1059 


0.8 


0.0947 


0.1429 


0.1174 


0.0918 


0.9 


0.0844 


0.1258 


0.1034 


0.0810 


1.0 


0.0762 


0.1118 


0.0922 


0.0726 






Table VI — 


3 = 8 




t 


Pcc(t) 


B + f . a (t) 


B + fi, a (t) 


Slit) 





1.0000 


1.0000 


1.0000 


1.0000 


0.1 


0.5756 


0.6335 


0.6088 


0.5842 


0.2 


0.4379 


0.5005 


0.4725 


0.4445 


0.3 


0.3727 


0.4256 


0.4006 


0.3756 


0.4 


0.3347 


0.3770 


0.3561 


0.3352 


0.5 


0.3099 


0.3432 


0.3262 


0.3092 


0.6 


0.2927 


0.3187 


0.3050 


0.2913 


0.7 


0.2802 


0.3005 


0.2895 


0.2786 


0.8 


0.2708 


0.2866 


0.2779 


0.2692 


0.9 


0.2637 


0.2760 


0.2690 


0.2621 


1.0 


0.2581 


0.2677 


0.2622 


0.2567 



This example will be used to show the operation of the error estimate 
(112). Using the increment h = 0.1, (158) and (159) were used to obtain 
0[e at fo, a (t)l * 2 [e ttt /o.«(t)] and 0[e«'/i.«(OL 2 [e*'/i.«(*)] at t = 0.5. 
Equation (122) was used to estimate € , a (0, €i,a(0- The error in si(t) was 
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Table VII — a = 12 



Pc(t) 



B+fo, a (t) 



B + fl,g(t) 



Slit) 






1.0000 


1.0000 


1.0000 


1.0000 


0.1 


0.6245 


0.6679 


0.6493 


0.6307 


0.2 


0.5255 


0.5629 


0.5458 


0.5286 


0.3 


0.4834 


0.5097 


0.4969 


0.4840 


0.4 


0.4611 


0.4789 


0.4698 


0.4607 


0.5 


0.4479 


0.4598 


0.4535 


0.4472 


0.6 


0.4396 


0.4476 


0.4427 


0.4378 


0.7 


0.4342 


0.4396 


0.4366 


0.4336 


0.8 


0.4306 


0.4342 


0.4322 


0.4301 


0.9 


0.4282 


0.4306 


0.4292 


0.4278 


1.0 


0.4265 


0.4281 


0.4272 


0.4262 



estimated by 2ci,„(£) - e , a (t) in which the estimates for e , a (t), e ha (t) 
were used. The results obtained are given in Table VIII. 



Table VIII — Error estimates at t = 0.5 

Estimate Si — / 



a 



tO,a 



Estimate 



«l,o 



Estimate 



4 0.0752 

8 0.0333 

12 0.0120 



0.0714 
0.0317 
0.0113 



0.0383 0.0372 0.0014 0.0030 

0.0163 0.0160 -0.0008 0.0002 

0.0056 0.0001 -0.0007 -0.0111 



VIII. SOME APPLICATIONS OF THEOREM 12 

The generating function, which will be designated G(z, t), of Theorem 
12, namely, 



G(z,t)=-^—f (- L -) 

1-2 \1 - 2/ 



(210) 



may sometimes be used to obtain explicitly the form of /„(£)• The fol- 
lowing are some examples. 
For f{t) = cos t, one has f (t) = 1/(1 + t 2 ), hence 

«*«>-arSfcr (211) 

The generating function for the Chebyshev polynomials, T n {t), of first 
kind is 15 



l-tz 



1 - 2*2 + 2 2 „=o 



= E T n (t)z n 



hence 



Gfet) = |„ z " (1 + £2) "'" +I,/2T "(^™)- 

One now obtains 
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(212) 



(213) 



r / t \21-(n+l)/2 

f n a)= L n cos t = yi + (^y) I T n 



V-fc) 2 . 



For /(i) = sin t, one has f (t) = t/(l + t 2 ), hence 
G(z, t) = l 



(214) 



(215) 



(l-z) 2 + t 2 

The generating function for the Chebyshev polynomials, U n (t), of second 
kind is 



l-2tz + z 2 n = 



= £ f/ n (t)z" 



(216) 



hence 



Thus 



G(z, t) = ± 2 n (l + t 2 )-("+»'HU n ( - 1 \ . (217a) 



/ n (t) = L n sin t = -4t [l + M-)* I 
n + 1 L \n + 1/ J 



-<n+l)/2 



xu n 



V-(^t) 2 . 



(217b) 



The Bessel functions provide additional interesting relations with 
orthogonal polynomials. For f(t) = Jo(t), one has/ (£) = l/Vl + t 2 , 
and 



G(z,t) = 



V(l-z) 2 + t 2 
The Legendre polynomials, P n (t), are generated by 

1 



Vl-2tZ + Z 2 n = 



= L P n (t)z n , 



(218) 



(219) 



hence 



1 + (— ) J Pn 



V-(^t) 2 . 

(220) 
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By the substitution of it for t, one derives immediately 

T / t \2"|-(n+l)/2 



P. 



V'-fcTT)' 



(221) 



Since Io(t) is convex, one also has 

I (t) < L n I {t) (222) 

for sufficiently large n. 
As another example relating to Bessel functions, consider f(t) = 

J Q {2Vt), then/ (t) = e~ l and 

G(2jt) = -J— c -t/(i-*>. (223) 

1-2 

The generating function for the Laguerre polynomials, L n (t), is 

_l_ e - t2 /(i-z) = £ L n (t)z» (224) 

1-2' n=0 

hence 

f n (t) = L„J (2v^) = e-*/<»+«L„ (jp^) • (225) 

IX. SUMMARY 

The methods of this paper have been found particularly useful in 
analyzing complex queueing phenomena whose Laplace transform 
representations are quite often implicitly defined. The error estimate 
of (112) has been found especially useful. Its computation is numerically 
effected by use of (158) and (159). 

It would be desirable to have an effective method of estimating the 
a parameter of (114) directly from /„ (t). In fact a method of this type 
which yields a rough evaluation has been devised and will be reported 
in a later paper. Of interest also would be further elaboration of the way 
structural properties of fit) are reflected in /„,«(£). 

The investigation of linear combinations of iterates, L r n , of the oper- 
ators L n may prove useful in providing additional enhancement methods. 
Especially, further investigation is needed concerning enhancement 
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methods which preserve the positivity of the approximation process. 

The isolated result of Theorem 16, which shows that exp(L n £ n f(t)) 
is a better approximation to f(t) than f n (t) when f(t) is log-convex, 
should be examined with the purpose of the possible construction of 
nonlinear approximation methods exploiting this structural charac- 
teristic. 
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APPENDIX 
Operations 

fit) 



fo(t) or f n (t) 



-fit) 

t 

f(at), a > 
e at f(t) 

Cf(x)dx 

fit) 
fit) 



tfit) 
tf{t) 



1 C't ( \ dx 

- I fo(x) — 
t Jo X 

fn(at) 

V " n+l) ' n \l-at/(n+l)) 

-J-.l f .( t L±±\ 

n+ly-o J \ n+l/ 

fo(t)-f(0) n + l\ in \1 

fo(t) - f(0) - tf(0) fi(t) - 2f (t/2) + f(0) 
t 2 ' t 2 



>D 



fnit)-2f n ^(^t) + f n .J r ^-t) 

(n + 1)2 Wl / U + 1 \n 



t 2 



i r l 

~ ( f(x)dx 
t Jo 

f(t)*h(t)(Me\\in) 



tfo(t) + t%(t) 

tfn(t) 

i r t 

7 I fn(x)dx 

t Jo 



fn(t)*h(t) 
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